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ABSTRACT 

Vortices have recently received much attention in the research of planet formation, as 
they are believed to play a role in the formation of km-sized planetesimals by collecting 
dust particles in their centres. However, vortex dynamics is commonly studied in non- 
self-gravitating discs. The main goal here is to examine the effects of disc self-gravity 
on vortex dynamics. For this purpose, we employ the 2D shearing sheet approxima- 
tion and numerically solve the basic hydrodynamic equations together with Poisson's 
equation to take care of disc self-gravity. A simple cooling law with a constant cooling 
time is adopted, such that the disc settles down into a quasi- steady gravitoturbulent 
state. In this state, vortices appear as transient structures undergoing recurring phases 
of formation, growth to sizes comparable to a local Jeans scale and eventual shearing 
and destruction due to the combined effects of self-gravity (gravitational instability) 
and background Keplerian shear. Each phase typically lasts about 2 orbital periods 
or less. As a result, in self-gravitating discs the overall dynamical picture of vortex 
evolution is irregular consisting of many transient vortices at different evolutionary 
stages and, therefore, with various sizes up to the local Jeans scale. By contrast, in 
the non-self-gravitating case, long-lived vortex structures persist for hundreds of orbits 
via merging of smaller vortices into larger ones until eventually their size reaches the 
disc scale height. Vortices generate density waves during evolution, which turn into 
shocks. This phenomenon of wave generation by vortices is an inevitable consequence 
of the differential character (shear) of disc rotation. Therefore, the dynamics of den- 
sity waves and vortices are coupled implying that, in general, one should consider 
both vortex and spiral density wave modes in order to get a proper understanding of 
self-gravitating disc dynamics. 

Our results suggest that given such an irregular and rapidly varying character of 
vortex evolution in self-gravitating discs, it may be difficult for such vortices to effec- 
tively trap dust particles in their centres, a proposed mechanism for planetesimal for- 
mation. Further study of the behaviour of dust particles embedded in a self-gravitating 
gaseous disc is, however, required to strengthen this conclusion. 

Key words: accretion, accretion discs - gravitation - hydrodynamics - instabilities 
- (stars:) planetary systems: protoplanetary discs - turbulence 



1 INTRODUCTION 

Research in protoplanetary disc dynamics is mainly fo- 
cused on investigating the dynamical activities of two basic 
types/modes of perturbations - spiral density waves and vor- 
tices. They provide a means of outward angular momentum 
transport necessary for the secular evolution of a neutral 
disc, where the magnetorotational instability cannot oper- 
ate. However, each mode is analysed for different settings: 
spiral density waves are the central perturbation types and, 
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therefore, the main agents responsible for angular momen- 
tum transport in self- gravitating discs, while vortices are 
commonly studied in non- self- gravitating discs. The forma- 
tion and dynamics of vortices in protoplanetary discs has 
only recently attracted attention because of the role they 
may play in planet formation. The latest developments have 
revealed that in non-self-gravitating discs vortices can be 
(linearly) coupled with spiral density waves. This notewor- 
thy finding, in turn, calls for revisiting self- gravitating disc 
dynamics with a particular emphasis on the possible role of 
vortical perturbations in the overall dynamical picture to- 
gether with spiral density waves. Below we first outline the 
recent results on the vortex dynamics in astrophysical discs. 
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It has been demonstrated in both 2D 
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of non- self- gravitating discs that only anticyclonic vortices 
(rotating in the same sense as the background shear) can 
survive in discs, though due to compressibility and/or 
viscosity they slowly decay on the timescale of several 
hundred orbital periods. Nevertheless, anticyclonic vortices 
can be thought of as long-lived structures. On the other 
hand, cyclonic vortices get strongly sheared by differential 
rotation of the disc and eventually disappear. Simulations 
are often initiated either with random potential vorticity 
(PV) perturbations, which contain positive (cyclonic) and 
negative (anticyclonic) values in equal portions, or with an 
imposed single vortex. The former shows that initial small 
scale regions where PV is negative (anticyclonic regions) 
first get sheared into strips and then start to wrap up into 
small scale vortices due to the specific instability discussed 
below. The latter gradually grow in size via merging 
into larger vortices until event ually their size becomes o f 
order the disc scale height (|johnson &: Gammid l2005h : 
growth beyond this scale is restricted by compressibility 
effects. The initial small scale cyclonic regions instead of 
wrapping up into distinct vortices get strongly sheared into 
strips and remain so during the entire course of evolution. 
This limitation on the vortex size was also confirmed in 
detailed siraulatio ns of a single vortex in a Keplerian disc 
(jBodo et al.ll2QQ7l ). In these simulations, a vortex with an 
initial length exceeding the disc scale height undergoes 
nonlinear adjustment, during which it decreases in size 
radiating excess energy in the form of spiral density waves 
and shocks, and finally settles down to scales equal to a 
few disc scale height. These final scales are independent of 
the initial vortex size and are determined by disc properties 
(sound speed and disc scale height). 

The emission of spiral density waves by vortices dur- 
ing the adjustment process is due to the background Ke- 
plerian shear that ensures (linear) coupling between these 
two modes when the horizontal scale of vortic es is equal to 
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or la r ger than the disc sca le height , 

2OO2I: I Johnson Gammid [2005 : Bodo et al.1 [2005. 2007; 

amatsashvili Chagelishvilj l2007l l . On the other hand, 
this condition is nothing more than the validity of 2D treat- 
ment. So, in the 2D case compressibility is a key feature 
that must be taken account of. In this respect, all the re- 
sults of 2D incompressible simulations of vortices should be 
regarded as approximate and less realistic. Emitted spiral 
density waves steepen into shocks afterwards and, hence, 
vortices appear to generate shocks that may play an impor- 
tant role in the outward transport of angular momentum, 
particularly in neutral discs wher e the magnetorotational in- 
stability cannot operate ( see e.g. iLi et al.ll200ll : lDavis|[200^ : 
[Johnson Gammid [iooa l. 

The basic underlying mechanism/instability responsi- 
ble for the development of vortices as w ell as the neces- 
sary criterion for that were identified by [Lithwickl (|2007[ l 
for incompressible shear flows. He interprets it as nonlinear 
Kelvin-Helmholtz instability of a vortex strip with vortic- 



ity of the same sign as the background vorticity (for incom- 
pressible flows vorticity and PV are equivalent); for a vortex 
strip with the opposite sign of vorticity no such instability 
is observed. This readily explains why cyclonic vortices can- 
not survive. This mechanism of vortex formation should be 
extended to compressible shear flows as well. Note that in 
this scenario one needs PV perturbations (e.g., vortex strips 
or random type) initially present in a disc. Other mecha- 
nisms of vortex, or PV, generation that do not necessarily re- 
quire this include inhomogeneities of entropy (temperature) 
dist ribution of a disc and are k nown as baroclinic instabil- 
itv ([ Klahr Bodenheime3[2003[ : iKlahr 2004; Peter sen et alJ 
2007] ) a nd R ossby wave instability (^Lovelace et alT " 19991 : 
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2006; Kla hr &: Bodenhe imer 20061). It has been shown that 
a smooth, sufficiently long-lived vortex is indeed able to ef- 
fectively trap dust particles in its core, possibly accelerating 
planetesimal formation. 

All the above-mentioned studies on vortex dynamics 
miss out an important aspe ct of protop laneta ry discs - self- 
gravi ty (s ee reviews by lAdams &: LinI Il993l . iDurisen et al.[ 
[2007[ and [Lodatd [2007[ on the role of self-gravity in pro- 
toplanetary disc dynamics). Typically in discs effective 
cooling times are too long to cause fragmentation unde r 
the action of their own self-gravity (Bolev et al.] [2006[ ). 
As a consequence, balance is established between heating 
due to gravi tational instability and cooling (self-regulation 
mechanism, [Bertin Lodatd [200 l[ l. Discs are expected to 
stay in this self-regulated quasi-steady gravitoturbulent 
state for a long time. In this case Toomre parameter 
hovers on the margin of gravitational instability. The 
properties of this state were nu merically investigated in a 
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Laughlin, Korchagin & Adams 
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20041, [2005[: [Mei ia et al.[ [2005[ : [Bolev et al 
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iStamatellos Whi tworth 20081, and references therein). A 
general dynamical picture is that spiral structure develops 
in a disc and transports angular momentum outwards 
through gravitational and hydrodynamic stresses, thereby 
allowing matter to accrete onto the central star. The 
angular momentum transport in this case is attributed to 
spiral density waves, which are thought of as the only per- 
turbation type present in the disc. In other words, almost 
all studies on self- gravitating disc dynamics concentrate on 
the dynamical activity of spiral density waves and leave 
another class of perturbations - vortices - out of consid- 
eration. As discussed above, the latter plays an important 
role in non-self-gravitating discs and it seems natural to 
look for them and analyse their nonlinear development in 
self- gravitating discs too. Indeed, in the linear theory it 
has already been shown that the coupling between spiral 
density wave and vortex m odes is even more efficient in 
the p resence of self-gravity (iMamatsashvili Chagelishvilj 
I2OO3). In perspective, this study will allow us to see if the 
same mechanism of planetesimal formation - dust particle 
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trapping by vortices - can also be at work in self- gravitating 
discs. 

In order to study the dynamics of vortical perturbations 
in self- gravitating discs, one must examine the behaviour of 
the PV field - the basic quantity characterizing vortex de- 
velopment. To date no systematic investigation of PV be- 
haviour, similar to that done for non- self- gravitating discs, 
has been carried out for self- gravitating ones. However, we 
should me ntion two relevant works b v Adams Watkins. 
(|l995l) andfWada. Meurer Normanl (2OO2 



Adams Watkinsl ([l995') investigated the dynamics of 



a single vortex in the quasi-geostrophic and local shearing 
sheet approximations in a self- gravitating Keplerian disc. 
The quasi-geostrophic approximation implies that the char- 
acteristic timescale of a problem is much larger than the 
orbital period. The vortex considered in their paper is in 
geostrophic balance and, therefore, remains steady for many 
rotation periods. The gas motion inside the vortex is sub- 
sonic as well. In this case the effect of self-gravity is only to 
make the effective length scale (Rossby radius) of the vor- 
tex larger than that in the non- self- gravitating case. The 
quasi-geostrophic approximation does not permit consider- 
ation of the most important aspect of dynamics - compress- 
ibility effects (spiral density waves and shocks), which are 
intertwined with vortices and have typical timescales of the 
order of orbital (shear) time. In other words, these relatively 
fast motions associated with compressibility are filtered out 
in t his approximation. 

IWada. Meurer NormanI (|2002l l investigated the prop- 
erties of the gravitoturbulent state in the interstellar 
medium of galaxies. The vorticity (but not PV) field in this 
state is indeed calculated in their paper, which has a rather 
complicated structure. However, the authors do not discuss 
in detail the properties of the vorticity field and its relation 
to the density and pressure fields. They only point out that 
negative and positive vorticity regions are associated with 
dense filaments seen in the density field. In that paper the 
main emphasis is placed on analysing the spectral proper- 
ties of gravitoturbulence. Their study clearly demonstrates 
that vortical perturbations are as important as spiral density 
waves in the formation of spectra of the resulting gravito- 
turbulent state. 

In the p r esent paper, following the appro a ch ad opted 
by iGammie (2001) and I Johnson Gammid (|2005h . we 
study the specific properties of vortex, or PV, evolu- 
tion in self- gravitating discs by means of numerical sim- 
ulati ons. We work in the 2D shearin g sheet approxima- 
tion ([Goldreich Lvnden-Bel l 1965; Go ldreich Tremaind 
1 19781 ) without invoking the quasi-geostrophic approxi- 
mation, thereby allowing for compressibility effects. In 
this respect our analysi s is more general than that of 
lAd ams Watkinsl (|l995l ). In addition, we do not impose 
a single vortex in quasi-equilibrium balance in the begin- 
ning. Instead, we start with random perturbations of veloc- 
ity components and, hence of PV, and trace the develop- 
ment of structures out of this chaotic field. These chaotic 
initial perturbations are more realistic than a single vor- 
tex. The main focus here is on the dynamical picture of 
PV evolution in the state of quasi-steady gravitoturbulence 
and how this picture differs from that occurring in non-self- 
gravitating discs given the same chaotic type of initial con- 
ditions. This, in turn, allows us to draw important conclu- 



sions about the effects of self-gravity on the formation and 
evolution of vortices. We will see that in the presence of 
self-gravity, vortex evolution is not as smooth and regular 
as it is in the non- self- gravitating case. The present study 
is actually an extension to the nonlinear regime with added 
cooling of iMamatsashvili Chagelishvilil ( 2007f ), where we 
carried out linear analysis of the transient growth and cou- 
pling of vortices and spiral density waves in self- gravitating 
discs. 

The paper is organised as follows. Physical approxima- 
tions and the mathematical formalism of the problem as well 
as the numerical techniques used are introduced in Section 
2, nonlinear evolution of the fiducial model and the com- 
parative study of correlation functions for different models 
are described in Section 3 and summary and discussions are 
given in Section 4. 



2 PHYSICAL MODEL AND EQUATIONS 

In order to study the evolution of vortices in thin self- 
gravitating gaseous discs, we employ the 2D local shearing 
sheet model. This model represents a simple analogue to a 
differentially rotating disc and has an advantage over global 
disc models in that it permits higher resolution study of dy- 
namical processes in discs. In the shearing sheet model, disc 
dynamics is studied in the local Cartesian coordinate frame 
corotating with the angular velocity of disc rotation at some 
radius from the central star. In this coordinate frame, the 
unperturbed differential rotation of the disc manifests it- 
self as a parallel azimuthal flow uo with a constant velocity 
shear in the radial direction. The unperturbed background 
surface density So and two-dimensional pressure Pq corre- 
sponding to this shear flow are assumed to be spatially con- 
stant. A Coriolis force is also included to take into account 
the effects of the coordinate frame rotation. As a result, in 
this local approximation the c ontinuity equation and equa- 
tions of motio n take the form (|Goldreich Tremaind [l 9 78l : 
iGammidlioOll ): 



- + V.(Eu)-,Qx- 



0, 



(1) 



This set of equations is supplemented with Poisson's equa- 
tion for a razor-thin disc 



AV^ = 47rG'(S-So)^(z). 



(4) 



Here u(ux,Uy), P,Y1 and are, respectively, the perturbed 
velocity relative to the background parallel shear flow 
uo(0, —qQx), the two-dimensional pressure, the surface den- 
sity and the gravitational potential of the gas sheet. Since 
equations (2-3) are written for perturbed velocities, only the 
gravitational potential due to the perturbed surface density 
S — So is used. In the dynamical equations, the gradients 
of the gravitational potential are taken at z = 0, i.e., where 
the shearing sheet is located, because only this quantity de- 
pends on the vertical coordinate z. Q is the angular velocity 
of the coordinate frame rotation as a whole, x and y are, re- 
spectively, the radial and azimuthal coordinates. The shear 
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parameter ^ = 1.5 for the Keplerian rotation considered in 
this paper. 

The equation of state is 

P = (7 - 1)U, 

where U and 7 are the two-dimensional internal energy den- 
sity and adiabatic index, respectively. We will adopt 7 = 2 
throughout. 

The central quantity of this study is the vertical com- 
ponent of potential vorticity referred to as PV for short: 



z-Vxu+(2 — 1 I duy dux 

dx dy 



+ {2-q)Q 



where z is the unit vector in the vertical direction. In the 
unperturbed state, where there is only the background Ke- 
plerian shear flow, Uo, with the constant equilibrium sur- 
face density So, PV is equal to /o = (2 - q)n/^o. The PV 
will play an important role in the subsequent analysis, as 
it generally characterizes the for mation of coherent stru c- 
tures (vortices) in a disc flow (e.^.,lGodon k Liviolll999al lbl: 
Ijohnson k Gammid [2005; "Bodo et al." l2007f ). Using equa- 
tions (1-3), after some algebra, one can show that the evo- 
lution of PV is governed by the following equation 



q^lx—- 
dy 



j_J_fd^dP_ as dP 
I dx dy dy dx 



This equation describes the advection of PV along the tra- 
jectories of fluid elements (Lagrangean derivative inside the 
brackets on the left hand side) and its change due to the 
nonlinear baroclinic term on the right hand side. In the 
present case, the pressure and surface density are not re- 
lated by any (e.g., isentropic, isothermal or polytropic) 
constraint, so this baroclinic term is, in general, non-zero 
and, therefore, PV is not conserved. However, in the lin- 
ear approximation it vanishes and PV is conserved making 
it possible to classify mo des into vo rtical and wave types 
([Mamatsashvih Chagelishvilil 120071 ). Also note that the 
gravitational potential, as it should be, does not explicitly 
enter into this equation; self-gravity only influences PV evo- 
lution through surface density, pressure and velocity flelds. 

The evolution of internal energy density is governed by 
the equation 

^ + V • (Uu) - qQx^ = -PV • u - - , (5) 

dt dy Tc 

where the flrst term on the right hand side is the compres- 
sional heating term and the second term takes into account 
cooUng of the disc. Here we assume the cooling time Tc to 
be constant and choose its value so that the disc does not 
fragment and achieves a saturated state, where all quantities 
fluctuate around constant average values. Namely, we take 
Tc = 20Q~\ which means a non-fragmenting disc according 
to Gammie's (2001) criterion. We refer to this quasi-steady 
state as gravitoturbulence. In the present study we concen- 
trate on examining the specific properties of PV evolution 
in such a gravitoturbulent state. Where fragmentation con- 
ditions are concerned, more realistic cooling laws, with Tc 
being a func tion of E, ?7, Q rather than a constant, are nec - 
essary (e.g., I Johnson &; Gammid l2003l : iBolev et al.1 12OO6I ). 

^ We also checked that as long as a disc does not fragment, differ- 
ent values of the cooling time do not qualitatively change results. 



Most simulations of self- gravitating discs with such a real- 
istic cooling, however, indicate that in this case cooling is 
usually not sufficient to cause fragmentation over most of 
the disc except possibly for the outer region s, therefore, a 
quasi - steady state is more likely (see review bv lPurisen et al.l 
I2OO7I ). We do not include any artificial heating terms in the 
internal energy equation; heating is solely due to the com- 
pressional term and to shocks captured by means of artificial 
(von Neumann- Richtmyer) viscosity in our code. 



2.1 Numerical methods 

Our computational domain in the {x,y) plane is a rectangle 
—Lx/2 < X < Lx/2, —Ly/2 < y < Ly/2 of size Lx x Ly 
divided into Nx x Ny grid cells. The numerical resolution 
is therefore Ax x Ay = Lx/Nx x Ly/Ny. We will assume 
that Lx — Ly = L and Nx — Ny = N. For the fiducial 
model presented below we take N = 1024, though we also 
ran a lower resolution (N — 512) model that converges 
to the fiducial model. In order to integrate the governing 
equations (1- 5) within this domain, we use a version of the 
ZEUS code (IStone k Normanlll992l). which is more suited 
to the shea ring sheet ([Gammid I2OOII : Ijohnson k Gammi3 
I2OOI I2OO5I ). ZEUS evolves these equations on a staggered 
mesh in a time-explicit, operator split, finite-difference fash- 
ion. As mentioned above, the code uses artificial viscosity to 
capture shocks. 

The transport scheme differs from that of the basic 
ZEUS algorithm. To show this, in the left hand side of equa- 
tions (1-5) we have explicitly decomposed the advection of 
physical quantities by the fiow into two parts: the first one 
(with u • V) represents advection by the perturbed veloc- 
ity and is done using the van Leer algorithm for calculating 
fiuxes of conserved quantities as in the original ZEUS al- 
gorithm. The second one (with —qQx • d/dy) describes the 
advection by the background Keplerian shear fiow up and is 
done by means of the FARGO scheme (jMassetl [ioooh with 
the use of the same van Leer algorithm for performing the 
fractional part of the shift; the integral part of the shift is 
done by a simple remap. This scheme has the advantage 
that it allows larger than standard integration time steps by 
removing the background shear flow from the Courant con- 
dition, which severely limits the time step because of large 
shear velocities at the x-boundaries. Now only the perturbed 
velocity u, rather than the full velocity uq + u, enters this 
condition. Numerical diffusion of the code is also reduced by 
this procedure. 

As with most works employing shearing sheet approx- 
imation, we adopt periodic boundary conditions in the y- 
direction and shearing periodic in the x-direction, that is, 
the x-boundaries are initially periodic but as time goes by 
they shear with respect to each other becoming again pe- 
riodic when tn = nLy/(qQLx) = n/(qQ), with n = 1,2,.... 
So, for each variable we can write: 

/(x, y, t) = /(x, y + L,t) {y boundary) 

f(x,y,t)=f{x-\-L,y — qQLt,t) (x boundary), 

where / = {ux,Uy, P,^,U). The way of implementing these 
shearing shee t boundary conditions is described in detail in 
iHawlev et all (|l995l ). The shearing in the x-boundary con- 
ditions, or shift in the ^-direction by an amount —qQ.Lt in 
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the ^^ghost zones^^, is done here by means of the FARGO 
algorithm. 

We define the autocorrelation function for PV as 
Ri{x,y) J / SI{x\y)6I{x + x\y + y)dxdy\ 

L/x-L/y J 

where the integration is over the entire rectangular simula- 
tion domain and 51 — I — Iq. The autocorrelation func- 
tion characterizes emerged coherent structures in a flow; 
its length scale can be identified with the characteristic 
scale of such struc tures. An analogous function was used 
by GammiCj (|200lh to analyse density structures in order 
to establish locality o f ang ular momentum transport and by 
Ijohnson &: Gammid |2005') to characterize coherent vortices 
in the non- self- gravitating shearing sheet. 

Our method of so lving Poisson's equati on is analo- 
gous to that used by Ijohansen et al. (l2007f) and differs 
from the origin al method of [ Gammid ( 200lh in the follow- 
ing respect. In IGammi the surface density is first 
linearly interpolated onto a grid of shearing coordinates 
{x,y = y - q^x{t - tp), tp = {qQ)-^NINT{qQt)), where 
it is exactly periodic, and then fourier transformed with a 
standard FFT technique. After that the fourier transform 
of the gravitational potential is found from the Poisson's 
equation rewritten in wavenumber (kx^ky) plane 



ip{kx,ky,t)\z^o = 



27iG^{kx,ky,t) 



where ^l^{kx,ky,t)\z^o and 12{kx,ky,t) are the fourier trans- 
forms of the gravitational potential and surface density with 
the background constant value subtracted. The potential is 
then transformed back into the grid of shearing coordinates 
and finally linearly interpolated onto the grid of {x,y) co- 
ordinates. We instead fourier transform directly from the 
{x,y) plane to the {kx,ky) plane and vice versa without 
doing an intermediate transformation to the shearing co- 
ordinates. However, in this case we take into account the 
fact that the radial wavenumber kx of each spatial fourier 
harmonic is no longer constant, but changes with time as 
kx{t) = ^(0) + qQkyt in order to remain consistent with 
the shearing sheet boundary conditions discussed above. As 
a consequence, at each time t in the corresponding compu- 
tational domain in the (kx^ky) plane there are only fourier 
harmonics with wavenumbers ky — 2'Kny/L^ kx — 2nnx/L-\- 
qQt' {27rny / L) , where t' = mod{t,l/{qQ\ny\)) {mod{a,b) is 
modulus after dividing a by b) and integer numbers Ux^riy 
lie in the interval —N/2 < Ux^riy < N/2. So, when do- 
ing fourier (inverse fourier) transform we decompose (sum) 
into (over) these wavenumbers. We accordingly modified the 
standard FFT routine to perform this decomposition (sum- 
mation). In some sense it is equivalent to changing from the 
{x,y) plane to the {x^y') plane via fourier transformation, 
which is more exact than that done by linear interpolation. 
In order to make the gravitational force isotropic on small 
scales, we keep only harmonics with k < 7tN/(L\/2), where 
k — {kx + ky)^^"^ . We have checked using the linear shearing 
wave test that the present version of potential solver gives 
more accurate results than that used previously. 



See website Ihttp:/ /imp.mcmaster.ca/ ~colinm/ism/rotfft.html| 



This version of the ZEUS code was tested on lin- 
ear problems involving shearing (vortical and com pressible) 
waves in the (non)- self-gravitating shearing sheet (Ga mmid 
I2OOII : [johnson Ga mmie 2005). The code results are in an 
excellent agreement with those of the linear theory until 
the decreasing radial wavelength of shearing waves becomes 
comparable to the grid cell size. So, we do not give these 
basic tests here and refer the interested reader to the above 
papers. 



2.2 Total energy equation and a parameter 

In this subsection we derive a relationship between the total 
energy and hydrodynamic and gravitational stresses. This 
relationship is important as it clearly shows the role of shear 
in the energy exchange between the background Keplerian 
flow and perturbations. The total energy density of per- 
turbed (from Keplerian shear flow) gas motion is defined 
as 

S= isu2+;7 + i(S-So)^, 

where the first term is the kinetic energy density, the second 
term, as defined above, is the internal energy density and the 
third term is the gravitational energy density. Using basic 
equations (1-5) and the shearing sheet boundary conditions, 
the equation for the average total energy is: 



Uy + 



4nG dx dy I Tc ' 



where the angle brackets denote averages defined as (/) = 
J fdxdy/LxLy^ where the integration is done over the en- 
tire simulation domain. The first and second terms inside 
the angle brackets on the right hand side are Reynolds and 
gravitational stresses respectively. As is clear from this equa- 
tion, these two terms are responsible for the energy exchange 
between the background Keplerian flow and perturbations, 
which is in fact possible due to the differential nature (shear) 
of Keplerian rotation. In the absence of background shear 
{q = 0) the total energy of perturbations changes only due 
to cooling. Incidentally, the sum in the angle brackets is 
also proportion al to the angular momentum flux density 
(|Gammidl200li ). Since the cooling time is not short enough, 
the system, as mentioned above, settles into a quasi- steady 
gravitoturbulent state in which d{E) / dt = (similarly for 
every average variable). Therefore, at these times we can 
write: 

or in terms of the well-known a viscosity parameter 



a : 



givmg 



1 /v. 1 
(Sci) \ 47tG 



1 



dip dip , 

^^^^ 
ox oy 



(6) 



^7(7 - l)VLrc ' 

where Cs is the adiabatic sound speed defined as = 
7P/S = 7(7 — 1)[//S. This relation is very important, as 
it shows that if a quasi- steady state can be achieved, then 
angular momentum transport must be determined solely by 
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the cooling time for a given rotation rate and, therefore, 
should be constant with t ime. (The same formula for a was 
derived bv 'Gam miell200ll in a different way.) It is the local 
analog of the result that the surface brightness of the disc 
does not depend on the details of the angular momentum 
transport mechanism. Again, the basic premise here is the 
possibility of the existence of a quasi-steady gravitoturbu- 
lent state, which is able to feed itself, or be self-sustained, 
at the expense of background shear flow energy extracted 
by Reynolds and gravitational stresses. In the presence of 
disc self-gravity and cooling this state is easily reached be- 
cause of the self-regulation mechanism (|Bertin Lodatol 
l200lh . However, in thin two-dimensional models of non-self- 
gravitating neutral discs steady outward angular momentum 
transport is hard to achieve. In this case typically the disc 
is unstable to vortex formation and develops well-organised 
vortices, which although able to transport angular momen- 
tum outwards mostly via generating compressible motions 
(shocks), are still slowly decayin g on the time scale of sev- 
eral hundred orbital periods (e . g . I Johnson &: Gammiell2005l : 
IShen. Stone Gardinerl l2Q06h . Thus, the role of self- gravity 
and cooling is crucial for the maintenance of gravitoturbu- 
lence. 

to the main analysis, we intro- 
variables. As mentioned above, 
state the background surface 
Pq and internal energy Uo are 
In the unperturbed state the 
7P0/S0 = 7(7 - l)Uo/^o. We 
variables: t Qt,{x,y) 
) {kxCso/Q,kyCso/Q), 
computational domain size 

) (Ux/CsO,Uy/CsO,Cs/Cso),T^ 



Before proceeding 
duce non-dimensional 
in the unperturbed 
density So, pressure 
all spatially constant 
sound speed is c^q = 
switch to non-dimensional 
{xQ/cso,yQ/cso), (kx.ky) 
correspondingly the 
L LQ/cso, {ux,Uy,c. 
E/So,(P,[/,^) - 



ip/c^Q,! IT^o/Q. Note that Cso/^ is actually the scale 
height of the disc. Therefore, distances are normalised by 
the disc scale height. The local Mach number is defined as 
M — y^Ux -\- ul/cs. The Toomre paramete r, a measure o f 
the self- gravity of a disc, is Q = cQ/ttGE (|Toomrelll96l ). 
These non-dimensional variables are used from now on 
throughout the paper. In the simulations presented below, 
we start with Q — Qq — Cso^^/ttCSo = 1.5. 



2.3 Initial conditions 

The initial conditions consist of random Ux and Uy pertur- 
bations superimposed on the mean Keplerian shear flow. 
The surface density and internal energy are not perturbed 
initially and thus are spatially uniform with values 1 and 
1/7(7— 1) respectively. To generate these initial conditions, 
for each velocity component at each point in the {kx,ky) 
plane we create a Gaussian random field with standard de- 
viation, or amplitude of power spectrum, given by the Kol- 
mogorov power law \ux,y{kx, ky)\'^ ~ k~^^^, k = yk^+k^, 
in the range kmin < k < kmax- The limits of this range are 
kmin = 27r/L and kmax — {N/ng)kmin, where Ug is the num- 
ber of grid cells contained within the smallest wavelength 
27r/ kmax for which we choose rig = 16 when N — 1024. 
Outside this wavenumber range the power spectrum is set 
to zero. This random field is then transformed back into the 
real [x^y) plane. For the amplitude of velocity perturba- 
tions a — (u^)"^/^, we initially take a — 0.6. Figure 1 shows 
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Figure 1. Initial (at t = 0) distribution of PV in the real plane 
corresponding to a Kolmogorov spectrum in wavenumber plane 
for the fiducial model with L = 20. Both negative (blue and light 
blue) and positive (yellow and red) values of PV are present. 



these initial conditions in terms of PV in the real plane. 
From this figure we can see that in the real plane the initial 
PV field has a chaotic character with equal contributions 
from negative and positive values. Although we start with 
a Kolmogorov power spectrum, a general dynamical picture 
emerging in the quasi-steady state does not depend strongly 
on the initial conditions. These random velocity perturba- 
tions are meant to mimic the turbulent state in a disc result- 
ing from the collapse of a molecular cloud (|Godon Livid 
2000). 

Note that we do not impose any constraints on the ini- 
tial fields Ux and Uy^ such as the requirement of incompress - 
ibility (V • u — 0) adopted b y Johnson k, Gamnii3 (|2005h 
and Shen . Stone Gardiner! (|20 06') to pick out onlv vor- 
tical perturbations; so these two velocity components are 
not correlated. We found that the mere presence of neg- 
ative PV values in initial conditions is sufficient for the 
development of vortices, even though they may also con- 
tain some frac tion of wave per turbations. Linear analysis 
( TMamatsashvili &: Chagelishvi i[2007) shows that in shear 
flows vortical perturbations become divergent if they are 
initially non-divergent and wave perturbations become ro- 
tational (with non-zero velocity curl) due to the background 
velocity shear. But wave perturbations never acquire PV if it 
is not present initially (of course, here we assume that there 
are no baroclinic terms in the PV equation in the linear 
approximation, or equivalent ly unperturbed surface density 
and pressure are uniform). Thus, truly vortical perturba- 
tions should be characterized by non-zero PV even though 
they may develop divergence in the course of evolution. It is 
the PV that determines the vortex formation and it must be 
present initially. By choosing initial perturbed velocity field 
with non-zero PV, we ensure the presence of vortical pertur- 
bations with a possible mix of wave perturbations with zero 
PV, which do not affect the vortex formation process. In 
other wards, the class of initial conditions capable of caus- 
ing vortex formation is broader; they must be character- 
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Figure 2. Evolution of the average kinetic, internal and grav- 
itational energies. During the burst phase, they undergo rapid 
(swing) amplification until t ~ 2 — 2.5 and then remain on av- 
erage constant in the quasi-steady gravitoturbulent phase, which 
sets in at about t = The evolution of the same quantities for a 
lower resolution {N = 512) run are similar. 



ized only by non-zero negative PV and may not necessarily 
satisfy the incompressibility requirement, that is, they can 
comprise a fraction of wave perturbations as well. This al- 
leviates the problem of vorticity injection. Initial conditions 
resulting from the collapse of a molecular cloud into a disc, 
where PV is everywhere zero would be somewhat contrived 
and therefore less likely. 



3 NONLINEAR EVOLUTION 

In this section we present the nonlinear evolution of the 
fiducial model with size L = 20 starting with the initial 
conditions described in the previous section. We describe 
in detail the development of structures in the PV, surface 
density and pressure/internal energy fields. The size L = 20 
corresponds to a minimum wavenumber ky min = 0.314 in 
the computation domain, for which linear swing growth is 
largest as a function of ky ai Q = 1.5. Thus, we ensure 
that from the outset our computational domain comprises 
those scales for which self- gravity is important. In the pres- 
ence of both strong Keplerian shear and self-gravity the 
main mechanism responsible for the growth of initial ve- 
locity perturbations is swing amplification instead of pure 
Jeans instability. Swing amplification has a transient expo- 
nential nature because of the ^Mrift^^ of radial wavenumbers 
of spatial fourier harmonics of perturbations through un- 
stable regions in the (kx^ky) plane that are broug ht about 
by the combined action of shear and self-gravitv (iToomre 
19811 : iKim Ostriked [2OO1I : iMamatsashvili Chage lishvilil 
20071 ) . As a result, the average values of various quantities 
undergo rapid amplification over a relatively short time in- 
terval (Figs. 2 and 3; see also Gammi 3 l200ll . The evolution 
of the same quantities for a lower resolution N = 512 run 
are almost similar. The small difference between these two 
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Figure 3. Evolution of (Q) and minimum Qrnin (upper plot) and 
the average Mach number (M), maximum Mmax and (cs) (lower 
plot). {Q) initially rises due to strong shock heating in the burst 
phase, but then levels off at about 3 in the gravitoturbulent phase. 
Qmin remains constant and low 0.6-0.7 implying that there are 
always unstable regions associated with negative PV, or vortices. 
(M) initially rises but then settles down to smaller values 0.3, so 
that gas motion is on average subsonic. However, Mm ax remams 
of order unity showing that there is still some sonic motion. The 
average sound speed (cs) also reaches a constant value in the 
gravitoturbulent phase. 



runs may be attributed to the specificity of PV, surface den- 
sity aii^_^nternal_ener^^ different resolutions, see 
also 'Shen. Stone Gard iner"2006'l . We call this interval the 
burst phase. During swing amplification, the initial velocity 
perturbations induce strong surface density perturbations in 
the form of shocks with density contrast by a factor of about 
100. These shocks have trailing orientation, because swing 
amplification always tends to produce trailing structures. 
The gas motion is also supersonic at this stage. It can be 
seen from Fig. 3 that the average and maximum Mach num- 
bers reach largest values 1.1 and 7 — 8 respectively. Shock 
regions are the sources of intense heating of the gas together 
with compressional heating (the first term on the right hand 
side of equation 5) and therefore the internal energy un- 
dergoes jumps in the shocks. Eventually strong shocks heat 
the disc up, but cooling compensates for heating, so that at 
about t = 5 a quasi- steady gravitoturbulent state is reached, 
where the average kinetic, internal and gravitational energies 
as well as (Q), Qmin, the average Mach number (M), the 
maximum Mach number Mmax and (cs) fiuctuate around 
constant values (Figs. 2 and 3). At these times, the angu- 
lar momentum transport parameter a is given by equation 
(6) and is constant. Now (M) fiuctuates around 0.3, so that 
on average the gas motion is subsonic, however, Mmax fiuc- 
tuates around 1.24 meaning that the motion is still sonic 
mostly in the vicinity of shocks. (Q) fiuctuates around 3, 
but we will see below that the Q{x,y,t) field is very inho- 
mogeneous containing values as small as 0.6-0.7 associated 
with negative PV regions (vortices). 

In Fig. 4 we trace the parallel evolution of PV, surface 
density, pressure (the same as internal energy for 7 = 2) 
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Figure 4. Evolution of PV, surface density, pressure/internal energy and Q parameter (from upper to lower rows respectively) for the 
fiducial model in the gravitoturbulent state. Coordinate axes are same as in Fig. 1. Snapshots are at times t = 14.7, 29.3, 44. The PV field 
consists of vortices with different sizes and strengths (blue and light blue regions). Vortices with sizes smaller than the local Jeans scale 
have higher (by absolute value) negative PV centres (blue dots). They also emit density waves during evolution that turn into shocks. In 
the surface density field, these vortices correspond to underdense central regions surrounded by overdense regions marking density wave 
emission places that are characterized by smaller (by absolute value) than the central values of PV. In the Q field, the underdense and 
overdense regions give rise to high and low Q values respectively, but the infiuence of self-gravity on such small scale vortices is weak. 
The other type of vortices are somewhat larger having scales comparable to the local Jeans scale. They are more diffuse with no clear 
shape (light blue regions in PV field) and are characterized by smaller (by absolute value) negative PV and stronger overdense regions 
(yellow and red in the surface density field) with even smaller Q than those in the previous case. For this reason, such vortices are in 
the process of being sheared and destroyed. In the pressure field, shocks and only the overpressure regions corresponding to stronger 
overdensities are noticeable fvellow and redV Exact values of these four quantities can be inferred from colour bars. 
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and Q fields to identify correlations between them and to 
compare them with those in non-self-gravitating discs. We 
emphasize from the beginning that the evolution is not as 
smooth and regular as it is in the non-self- gravitating case. 
As we will see below, in self- gravitating discs vortices are 
transient structures undergoing recurring phases of forma- 
tion, growth and destruction, whereas in non- self- gravitating 
discs small scale vortices once formed, gradually merge into 
larger vortices. 

During the burst phase, initially small scale positive 
and negative PV regions get strongly sheared into strips, 
but negative PV (anticyclonic) regions start to wrap up into 
vortex-like structures towards the end of the burst phase du e 
to nonlinear Kelvin-Helmholtz instability ([Lithwickl l2007l ) . 
The positive PV regions remain sheared into strips showing 
no signs of vortex formation during the entire course of evo- 
lution. Thus, only anticyclonic regions are able to survive 
in shear flows by taking the form of vortices. As mentioned 
in the Introduction, this is a quite general result confirmed 
in many simulations of astrophysical as well as purely hy- 
drodynamical contexts. As evident from Fig. 4, in the sub- 
sequent quasi-steady state the overall dynamical picture of 
PV evolution is still irregular and varying on the dynamical 
timescale. For this reason, throughout the paper we use the 
term ^Vortex^^ in a broader sense meaning generally a neg- 
ative PV region and not only a vortex with a well-defined 
shape commonly occurring in non-self-gravitating discs. As 
demonstrated below, in a quasi-steady state these negative 
PV regions are associated with underdense/underpressure 
and overdense/overpressure regions superimposed on the 
shocks that may be important for trapping dust particles. 

Small scale anticyclonic vortices initiated at the end of 
the burst phase continue to grow further in size in the quasi- 
steady saturated state which is reached, as mentioned, at 
about t = 5. Figure 4 presents a typical dynamical picture 
of such a gravitoturbulent state, which is quasi- stationary 
and remains unchanged (on average) with time during the 
entire course of evolution. Vortices with a large (by abso- 
lute value) negative central value of PV (blue dots in the 
PV field) have more or less vortex-like shape and corre- 
spond to the central underdense regions surrounded by the 
higher density regions. Due to the background Keplerian 
shear, the total motion in the vicinity of such negative PV 
regions is a complex mixture of vortical and compressive 
(wave) motions giving rise to this type of density structure. 
In these enhanced density regions, PV is smaller by abso- 
lute value than that in the vortex centre. Thus, generally in 
shear fiows the coupling with wave motions is typical of vor- 
tex dynamics irrespective of the disc being self-gravitating 
or not. These compressive motions turn into shocks after- 
wards that are also seen in the surface density and pres- 
sure fields (more precisely, the shock structure in these fields 
is a net result of nonlinear superposition/mixture of shock 
waves emitted by such individual vortices). So, classifying 
perturbations into purely vortical - having only zero diver- 
gence and non-zero curl - and compressive - having only 
non-zero divergence and zero curl - is not quite appropri- 
ate, because due to background shear, divergence and curl 
become linked and change with time during the evolution 
that, in turn, causes mixing of these t wo perturbation types 
(iMamatsashvili Chagelishvilill2007l l. With regard to an- 
gular momentum transport, it would be more exact to at- 



tribute it to the overall gas motion with non-zero PV and 
divergence rather than to compressive (wave) or to vortical 
perturbations only. 

A similar dynamical picture is also observed in the evo- 
lution/adjustment of a single v ortex in a compre ssible non- 
self- gravitating Keplerian disc teodo et al.ll2007l l. The vor- 
tex initially produces saddle- like density structure, much like 
observed in our simulations, with an underdense region at 
the location of the vortex centre surrounded by an overdense 
region with accompanying shock waves . 

With time the effect of self-gravity on each vortex 
comes into play because of the growing length scale of vor- 
tices and at the same time of the increasing surface den- 
sity /overdensity associated with them. The vortex growth 
in size is a consequence of inverse - towards larger scales 
- cascade of energy characteristic of 2D flows as clearl y 
demonstrated in simulations (e.g.. [Codon Livid Il999al). 
The ampliflcation of the vortex surface density is again 
due to the swing mechanism. For the same reason that in 
self- gravitating discs perturbations experience considerably 
larger swing growth than in non-self-gravitating discs, the 
traces of vortices in the surface density field - the overdense 
regions around underdense regions - become more and more 
noticeable on the background of density variations related 
to shocks. From Fig. 4 it is seen that overdensities are also 
characterized by small Q values. In the pressure/internal en- 
ergy field we see shocks and only overpressure regions. The 
stronger overdensities (yellow and red in Fig. 4) associated 
with the final stage of vortex life (see below) clearly cor- 
relate with these overpressure regions. However, it is still 
hard to find as good a correspondence/similar features be- 
tween PV and pressure/internal energy fields as it is for the 
surface density field. Nevertheless, in self- gravitating discs 
one can identify distinct features in the surface density and 
pressure/internal energy fields connected with vortices. By 
contrast, in analogous non- self- gravitating isothermal simu- 
lations the situation is different. As mentioned above, similar 
density structures, i.e., with overdense and underdense re- 
gions, were also observed in numerical simulations of a single 
vortex before its settlement into an equilibrium configura- 
tion, where the a nticyclonic vortex gives rise only to a single 
overdense region (|Bodo et al.ll2007l ). However, other simula- 
tions of vortex formation from an i n itial random distribution 
of PV ( J ohnson Gammid l2005l : [Sh en. S tone Gardiner! 
2006), which are not limited to only a single vortex, did 
not find overdense regions easily identifiable with individual 
vortices in the surface density field; there were only varia- 
tions in the surface density due to the shocks shed by these 
vortices (see also Fig. 5). In our case instead, swing amplifi- 
cation due to self-gravity makes overdense and underdense 
regions clearly noticeable on the background of shocks. 

Since Q < 1 at the location of the stronger overdense 
regions, self- gravity becomes dominant in the dynamics of 
vortices and appears to strongly shear and deform the latter 
over about 1 — 3Q~^ and most importantly to inhibit their 
further growth in size. In other words, self-gravity opposes 
an inverse cascade of energy to smaller wavenumbers and 
scatters it to higher wavenumbers. This fact is also refiected 
in the autocorrelation function of PV shown below (Fig. 6). 
The vortices can only grow to a size comparable to the local 
Jeans scale Aj ^ Qcs (in non-dimensional variables). When 
the vortex size approaches the local Jeans scale, it has no 
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Figure 5. Typical PV field (upper plot) at f = 44 for the isother- 
mal fiducial model with self-gravity switched off. The sound speed 
is equal to the average sound speed in the gravitoturbulent state 
(see Fig. 3). Also shown is the corresponding surface density field 
(lower plot). It is hard to see traces of individual vortices in the 
surface density field; only shocks shed by these vortices are visible. 
Vortices are organised into well-shaped larger scale structures as 
opposed to what is observed in the self-gravitating case in Fig. 4. 



clear centre with high (by absolute value) negative PV re- 
gion surrounded by smaller (by absolute value) PV region 
and, therefore, corresponds only to the stronger overdense 
region (yellow and red regions in the surface density field in 
Fig. 4) with no underdense centre. At this time the shape 
of the vortex is more irregular, sheared and diffuse charac- 
terized by about 10 times smaller (by absolute value) PV 
than that of growing vortices. From the pressure/internal 
energy field it is clear that in these stronger overdense re- 
gions the internal energy is higher as well and it is expected 
that corresponding Q will eventually rise and switch off self- 
gravity /gravitational instability at that location. Then the 
overdense region will quickly disperse, or get sheared away. 
After that sheared PV regions with higher Q, in turn, can 
start another cycle of vortex formation by undergoing the 



same nonlinear Kelvin-Helmholtz instability and again wrap 
up into vortices. Thus, in self- gravitating discs vortices have 
recurring nature - each vortex forms, grows to a size of the 
order of the local Jeans scale and after that gets sheared and 
destroyed by the combined effects of self-gravity and Keple- 
rian shear. Each such a cycle typically lasts for about 2 or- 
bital periods or less (1 orbital period= 2ti /Q). Then sheared 
PV regions turn into vortices and after that the whole pro- 
cess starts again. The cooling ensures that the minimum 
Qmin is kept low (0.6 — 0.7) and nearly constant with time, 
so that self- gravity continues to play a role. 

Above we have described the behaviour of some indi- 
vidual vortices only. From Fig. 4 we see that the PV field 
contains vortices with different sizes that never get organ- 
ised into distinct coherent vortices as they do in the non-self- 
gravitating case (Fig. 5). Vortices evolve differently - some 
of them are at the end of evolution having already grown to 
the Jeans scale sizes and are characterized by clearly iden- 
tifiable overdense regions and low Q, but others have yet to 
go through this phase and, therefore, may still have under- 
dense central regions surrounded by higher density regions 
and correspond to higher Q. Thus, regions of high and low 
Q coexist throughout the disc at all times. A similar gravi- 
tot urbulent state was observed in th e global disc simulations 
by IWada. Meurer NormanI (l2002l ) in the context of ISM 
turbulence. However, these authors did not measure PV cor- 
responding to high and low Q regions. In general, the overall 
dynamical picture, as mentioned, is very irregular and it is 
hard to keep track of individual vortices for several orbital 
periods, because they are short-lived structures. Figure 5 
shows the stark contrast between non-self-gravitating and 
self- gravitating cases. In this figure we present a snapshot of 
the same fiducial model, but without self-gravity and with 
constant sound speed (isothermal equation of state) equal 
to the average sound speed of the fiducial model in the 
gravitoturbulent state (see Fig. 3). In this case the dynami- 
cal pi cture is identical to that described by other authors 
fe.g.. lUmurhan k Regeyl l2004l : I Johnson Gammid l2005l : 
IShen. Stone Gardinei]l2006h . Now it is easier to trace each 
vortex and to see how they merge into larger vortices. In 
the self- gravitating case, irregular and chaotic phases of vor- 
tex evolution is quasi-steady, or equivalently self-sustained 
during the entire course of evolution, whereas in the non- 
self- gravitating case, though vortices are well-organised and 
regular, they gradually decay. 

We also carried out simulations for two other models 
with sizes L — A and L = 10 and with the same type of 
initial conditions. The dynamical picture of vortex evolu- 
tion described in detail for the fiducial model remains qual- 
itatively the same. There are differences only in saturation 
times, and (average) values of various quantities in the grav- 
itoturbulent state. In Fig. 6 we compare the autocorrelation 
functions of PV for these models with that of the non-self- 
gravitating isothermal model plotted in Fig. 5. The extent 
of the autocorrelation function, or the correlation length, 
remains on average unchanged with time in the gravitotur- 
bulent state and does not depend on the spectrum of initial 
conditions. For the given model size L — 20, the correla- 
tion length is smaller than that in the non-self-gravitating 
case implying that vortices are less coherent in the gravito- 
turbulent state (we remind the reader that the sound speed 
of the isothermal model, determining the correlation length 



Vortices in self- gravitating gaseous discs 11 




2 10 12 -2-1012 




X X 



Figure 6. Autocorrelation functions of PV for four models: 
isothermal non-self-gravitating with L = 20 (Fig. 5) and three 
self-gravitating with L = 20, 10, 4 in the gravitoturbulent state. 
For the two fiducial models with the same L = 20 the correlation 
length is largest for the isothermal non-self-gravitating model and 
greatly decreases in the presence of self-gravity. It also decreases 
with decreasing L. 

for this model, is equal to the average sound speed of the 
fiducial self-gravitating model) . This again can be explained 
by the tendency of self-gravity to oppose inverse cascade 
of energy to larger scales by scattering it to smaller scales 
and thereby broadening the spectrum of the PV autocor- 
relation function. An alternative explanation is that in the 
self- gravitating case the contribution of Jeans scale size PV 
regions in the correlation function is smaller than that of 
smaller size PV regions. The correlation length appears to 
decrease with decreasing model size, whereas in the isother- 
mal case it does not depend on the si ze of the model and 
is det ermined by the disc scale height ( Johnson Gammid 
In any case, correlation length is always smaller than 
the model size justifying the local treatment of vortices. 



4 SUMMARY AND DISCUSSIONS 

We have studied the specific properties of potential vor- 
ticity (PV) evolution in self- gravitating discs in the shear- 
ing sheet approximation. The evolution has been traced via 
numerical integration of the basic hydrodynamic equations 
supplemented with Poisson's equation taking care of self- 
gravity. Since we are interested particularly in the proper- 
ties of vortex evolution in quasi-steady gravitoturbulence, 
we have chosen a simple cooling law with a constant cool- 
ing time large enough so that the disc settles down into this 
state. Our analysis has shown that in the self- gravitating 
case, vortices appear as transient structures undergoing re- 
curring phases of formation, growth to sizes comparable to a 
local Jeans scale, and eventual shearing and destruction due 
to the combined action of self-gravity (gravitational insta- 
bility) and background Keplerian shear. Each such a phase 
lasts for a few orbital periods. As a result, the overall dynam- 



ical picture is irregular, consisting of many transient vortices 
at different evolutionary stages and, therefore, with various 
sizes up to the local Jeans scale. By contrast, in the non-self- 
gravitating case, long-lived vortices form and grow in size via 
merging into larger ones until eventually their size reaches 
the disc scale height. The motion within vortices, or more 
precisely in the vicinity of negative PV regions, is a complex 
mixture of compressive (wave) and vortical motions which 
are difficult to separate from each other. Compressive mo- 
tions turn into shocks afterwards, or equivalently, vortices 
appear to generate shocks. Therefore, the dynamics of com- 
pressive motions (density waves) and vortices are coupled 
implying that, in general, one should consider both vortex 
and spiral density wave modes to get a proper understanding 
of self- gravitating disc dynamics. 

It is well known that overpressure/overdensity regions, 
if pre sent in a disc, can a ct as t raps for d ust particles 
(e.g.. iHaghighipour k [20031: iRice et al.1 [2004 l2006l : 
iFromang k NelsonI l2005l : iLvra et all l2008h . As has been 
demonstrated in our simulations, in self- gravitating discs 
anticyclonic vortices, or generally negative PV regions, are 
able to produce quite noticeable overpressure/overdensities, 
though not long-lived. They have a transient character and 
vary on the dynamical timescale. Hence, given such an irreg- 
ular and rapidly changing nature of vortex evolution in self- 
gravitating discs, it seems difficult for corresponding over- 
densities/overpressures to effectively trap dust particles in 
their centres. Further study of a coupled system - dust par- 
ticles embedded in a self- gravitating gaseous disc - is, how- 
ever, required to strengthen this conclusion. 

Here we have considered the simplified case of a razor- 
thin (2D) disc in order to gain first insight into the effects of 
self- gravity on vortex dynamics. Obviously, for a fuller un- 
derstanding 3D treatment is necessary. The situation can be 
different in the 3D self- gravitating cas e, which is evidently 
stratified in the vertical direction. iBarranco MarcusI 
(2005; ) demonstrated that in the 3D non- self- gravitating 
stratified shearing box, off-midplane vortices appear and 
survive for many orbita l perio ds. A similar 3D study by 
IShen, Stone Gardiner! (|2006h , but not including stratifi- 
cation, showed that vortices are unstable and get quickly 
destroyed due to an elliptical instability. So, it can be said 
that the presence of stratification helps the survival of vor- 
tices. On the other hand, as we have seen here, self-gravity 
does not favour long-lived vortices. But in the 3D case the ef- 
fect of self-gravity is somewhat reduced compared with that 
in the 2D case considered here. Thus, in a 3D generalization 
of the present problem there are two competing factors - 
stratification and self-gravity - and it remains to be seen in 
numerical experiments which of these two prevails. Another 
question of interest in the 3D c ase (either self-^avit a tin^ o r 
not) is, as pointed out also by Ijohnson Gammid (|2005[ ). 
if long-lived off-midplane vortices can be generated from a 
random PV distribution, as happens in the 2D case, rather 
th an from specially chosen vo rtex solutions as in simulations 
of iBarranco k MarcusI (|2005l ). 
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